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Abstract 

The bifurcation diagram of a truncation to six degrees of freedom of the equations for 
quasi-geostrophic, baroclinic flow is investigated. Period doubling cascades and Shil'iiikov 
bifurcations lead to chaos in this model. The low dimension of the chaotic attractor suggests 
the possibility to reduce the model to three degrees of freedom. In a physically comprehensible 
limit of the parameters this reduction is done explicitly. The bifurcation diagram of the reduced 
model in this limit is compared to the diagram of the six degrees of freedom model and agrees 
well. A numerical implementation of the graph transform is used to approximate the three 
dimensional invariant manifold away from the limit case. If the six dimensional model is 
reduced to a linearisation of the invariant manifold about the Hadley state, the Lorenz-84 
model is found. Its parameters can then be calculated from the physical parameters of the 
quasi-geostrophic model. Bifurcation diagrams at physical and traditional parameter values 
are compared and routes to chaos in the Lorenz-84 model are described. 

1 Introduction 

The equations for atmospheric flow form one of the most intensely studied dynamical systems of 
the last century. Both their practical importance and their mathematical richness have attracted 
much attention. The atmospheric equations are studied in various forms, depending on the physical 
domain, and the time and length scales of interest. The starting point of the analysis in this paper is 
a model which describes midlatitude atmospheric flow on a synoptic scale, a few thousand kilometers 
in space and a week or so in time. A phenomenological description of the typical flow patterns in 



this range can be found in Peixoto and Port [1992 1 , chapter 7. The model consists of the filtered 



equations, derived from the basic atmospheric equations under the assumption of quasi-geostrophic 



(QG) and hydrostatic balance (see Holton [1992Q 



In the absence of dissipative processes and forcing through solar heating, the filtered equations 
are energy preserving. The dominant dissipative terms are friction at the earths surface, internal 
friction and Newtonian cooling. The solar heating induces a strong temperature gradient in the 
meridional direction. Additional temperature gradients in the zonal direction can be induced by, 
e.g. land-sea contrast. The response to this forcing is a strong westerly circulation, called the jet 
stream. This circulation can become dynamically unstable so that traveling waves develop. 

The jet stream pattern is nearly equivalent barotropic, which means that its height dependence 
can be represented by multiplication by a scale function. In other words, at each surface of constant 
height, the velocity field has the same shape, but may have a different amplitude. In contrast, the 
traveling waves can be baroclinic, which means that their phase depends on the vertical coordinate. 
Typically, they exhibit a westward tilt with height. Theoretical studies of the filtered equations 
indicate, that the baroclinicity of these traveling waves changes in the course of their life cycle 
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[Frisius, 199?]. In the growing, strongly baroclinic, phase, they extract energy from the jet stream. 
In the decaying phase, they can become equivalent barotropic and transfer energy back into the jet 
stream. 

The focus of this study is on the interaction of the jet stream and the baroclinic waves and its 
representation in a low order model. With the aid of discretisation in the vertical and Galerkin 
truncation in the horizontal coordinates, we approximate the filtered equations by a finite number 
of ordinary differential equations (ODE's). The discretisation can be done with out violating the 
energy preserving nature of the filtered equations, as described in Lorenz [196C |. The number of 
layers is fixed to two, the minimal number necessary to describe baroclinic waves. 

The two layer model is considered on an /-plane, a rectangular domain on which the Coriolis 
force is taken to be constant. The partial differential equations of the two layer model are then 
projected onto Fourier modes. In each layer, we use one zonally symmetric pattern, representing 
the jet stream, and two patterns which combine to represent a traveling wave. Thus, the ODE 
model has six degrees of freedom. The solar forcing is represented by constant terms, which are 
used as bifurcation parameters. 

The bifurcation diagram is organised by its codimcnsion two points, namely fold-Hopf, 2:1 
resonance and a neutral saddle-focus on a homoclinic bifurcation line. Two routes to chaos are 
readily identified: period doubling cascades and a Shil'nikov type bifurcation. An inspection of the 
spectra of equilibria and periodic orbits found, and the calculation of the Kaplan- Yorke dimension of 
the chaotic attractor, leads to the conjecture, that there is a three dimensional, globally attracting, 
invariant manifold in the phase space of the six dimensional model. 

In order to calculate a first approximation of this invariant manifold, a small parameter is intro- 
duced into the equations. In the limit where this parameter tends to zero, an analytic expression is 
obtained. This limit has a clear physical interpretation in terms of the energy transfer between the 
jet stream pattern and the traveling waves. Numerical evidence for the persistence of the invariant 
manifold away from this limit is obtained by the use of techniques introduced by Broer et al. [1997 1 
and [Foias et al. [1988| ]. 

Reducing the six dimensional model to the invariant manifold, in the limit where the small 
parameter tends to zero, a three dimensional model is obtained. Its bifurcation diagram is compared 
to the corresponding diagram of the six dimensional model in order to see if the qualitative dynamics 
is retained. This comparison is convincing. Particularly, the codimension two points are still present. 

If the six dimensional model is reduced to a linearisation of the invarian t manifold ab out the 
Hadley state, the Lorenz-84 model emerges. This model was introduced by Lorenz [1984 1 as the 
simplest model capable of representing the basic features of midlatitude, synoptic flow. To the 
author's knowledge, no derivation of the Lorenz-84 model fr om atmospheric flo w equations has 
been presented before. A rather ad hoc link was established by Wiin-Niclscn [1992 1 , but in his work 
the reduction to three degrees of freedom is not based on physical or mathematical arguments. The 
link established here enables us to calculate the parameters in the Lorenz-84 model from the physical 
parameters in the filtered equations. As it turns out, one of the parameters comes out significantly 
different from its traditional, yet unmotivated, value. A continuation in this parameter relates the 
bifurcation diagram found at the traditional parameter value, presented in Bhilnikov et al. [1995| ] , to 
the one found at the physical value. The latter still bears resemblance to the bifurcation diagram 
of the six dimensional model, but the neutral saddle-focus transition is no longer there. Hence 
the route to chaos through a Shil'nikov type bifurcation is absent. It is shown, that chaos through 
period doubling cascades, the Ruelle-Takens scenario and intermittency does occur in the Lorenz-84 
model. 

The derivation presented here is not unlike the derivation of the Lorcnz-63 model from the fluid 
dynamical equations governing Rayleigh-Benard convection [Saltzman, 1962: Lorenz, 1963]. There 
too, the Galerkin truncation calculated has six degrees of freedom and can be reduced to a three 
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dimensional invariant manifold. In that case, however, the invariant manifold is linear. 

Galerkin truncations of the filtered equations have been studied on various domains and at 
various truncation numbers, see De Swart [1988 and references therein. One lesson to be learned 
from these studies, is that severe truncations, such as the one studied here, can only be regarded as 
qualitative models. A quantitative comparison to solutions of the filtered equations may be sensible 
at a trunc ation number in the order of a hu ndred or higher, depend ing on the basis functions used 
(see, e.g., Achatz and Branstator [1999 and Itoh and Kimoto [1996 1 ) . Low order models, however, 
allow us to isolate a physical process , such as the interaction between the jet stream and baroclinic 
waves, and represent it in a simple way. Because of this conceptual simplicity, and the fact that 
they are easy to integrate numerically, low dimensional truncations are widely used for testing and 
illustrating new ideas in dynamical systems theory, meteorology and climatology. The L orcnz- 
84 model for instance, has been used to investigate low- frequency atmospheric variability [Piclkc 



and Zcng, 1994 1, me asures of predic t ability | Gonzalez-Miranda, 1997 and time scale interaction in 



the climate system |Roebber, 1995; van Veen et al., 2001 1. The link with the filtered equations 



presented here, validates the use of the Lorenz-84 model in these contexts, though with different 
parameter values. 



2 The QG two layer model 



As a starting point for the calculations we take the quasi-geostrophic two layer model, described by 
Lorenz [1960|. Alternative derivations of this model can be found e.g. in Holton [1992| , chapter 8, 
or the review article by De Swart [1988|. In Lorenz' article the stress is on the energy conserving 



nature of the nonlinear interaction terms, in De Swarts review article strong scaling arguments are 
provided. 



2.1 Setup of the model 

In the quasi-geostrophic approximation, the dry atmosphere is described by the velocity field, v, 
and the temperature field, T. It is convenient to use the streamfunction, the velocity potential, 
X, and the potential temperature, 0, as variables. As a vertical coordinate we use pressure instead 
of height. The velocity and the temperature can then be expressed as 

v r =kx V* v d = Vx 

v = v r + v d + cjk T = e(j^j " P (1) 

Here v r and are the divergence free and the irrotational part of the horizontal velocity, p is 
pressure, p s is surface pressure, to = dp/dt is the vertical velocity, k is the vertical unit vector and 
c v and c p are the specific heat of dry air at constant volume and pressure, respectively. 

At the earth's surface, the lower boundary, we impose that p — p s is constant and w — 0. At the 
upper boundary we have p = and lu = 0. Discretisation of the vertical in layers means that we 
replace each function of three spatial variables by a number of functions of longitude and latitude 
only. In the simplest case we take two layers, the minimal number necessary to describe baroclinic 
waves. 

The two layers are bounded by three isobaric surfaces at po — p s , p2 = p s /2 and p4 = (see 
figure Q). Vertical derivatives are replaced by linear interpolations, e.g. in the continuity equation: 

V X+ dp~-° { V 2 X i-c(p 2 ) = ^V( X i+X3)-Q (2) 



3 



p 4 =o 



Layer 2 



p^= - p 

n 4 



P 2 = i P S 



Layer 1 



//////////////////// 



Po=P s 



Jgure 1: The vertical discretisation in pressure coordinates. If the static stability, <r, is fixed, equations (4) and 
3) determine t, and xi- 



The streamfunction and the potential temperature in the lower and the upper layer are denoted by 
$ 13 and Oi,3, respectively. The pressure in the lower and the upper layer is set to p\ = 3p s /4 and 
j>3 = p s /4. The equations are written in terms of vertical means and differences, defined as 

\& = l/2(5 , 3 + ^i) the barotropic streamfuction, 

r = 1/2(^3 — ^x) the baroclinic streamfuction, 

= 1/2(03 + ©l) the mean potential temperature, 

(7=1/2(63-61) the static stability. (3) 

The static stability, a, will be taken constant. 

In addition to the conservative dynamics described in Lorcnz [1960[ , we introduce linear damping 
through friction at the earth's surface, the terms proportional to C, friction at the boundary of the 
two layers, the term proportional to C and Newtonian cooling, the term proportional to hjq. The 
temperature forcing is given by 6*. The resulting equations are 



— V 2 * = -J(tf , V 2 * + /) - J(r, V 2 t) - CV 2 (* - t) (4.1) 

^-V 2 r = -J(r, V 2 * + /) - J(*, V 2 r) + V • (fVxi) + CV 2 (* - r) - 2C"V 2 r (4.2) 
at 

|U = -</(*, 6) + ( tV 2 xi-M©-©*) (4-3) 
ot 

where the Jacobian operator, J, is defined as J{A, B) — V A ■ V ' B x k for any pair of functions 
A, B. The Coriolis parameter has been denoted by /. 

Furthermore we have the thermal wind equation, relating the shear streamfunction to the mean 
potential temperature 

6c p V 2 6 = V • (/Vt) (5) 
Where b = \ [(p^j " P ~ Op) ] ~ 0.124 comes out of the discretisation scheme in the vertical. 
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2.2 Domain and boundary conditions 



The set of equations (4) will be considered on an /-plane, a rectangular domain centered about a 
fixed latitude <po on which the Coriolis parameter is approximated by the constant value fo. This 
domain has length L in the zonal direction and D in the meridional direction. On this plane we 
will use Cartesian coordinates x £ [0, L) and y £ [0, D]. In the following we set <f>o = 45°. 

In the zonal direction we take periodic boundary conditions. In the meridional direction we 
have 



<9<3> dr dxi,3 
dx dx dy 







at y = 0, D 



(6.1) 



This mean s that there i s no mass flux through the boundaries. The second condition was put 
forward by Philips [1954 1, and imposes that there is no net flow along the boundaries. 



o dy J dy 



at y = 0, D 



(6.2) 



It follows from the thermal wind equation (|5|), and the restriction that there be no net heat flux 
through the boundaries, that satisfies 



dx 







50 
dy 



dx = Q at y = 0, D 



(6.3) 



With these boundary conditions we can consider to describe deviations from the spatially 
averaged potential temperature. 

2.3 Scaling 

In table ([j]) the scales, suitable for the synoptic physics are listed. In the right column the numerical 
values, used below, are listed. In the following all quantities are dimensionless, unless otherwise 
indicated. The dimensionless length of the domain will be denoted by s = L/D. 



Length 


D 


5 • 10 3 km 


Time 




7 days 


Temperature 


n - bc p 


34.3K 


Mass 


M = p f 6 R 


7.7 ■ 10 18 kg 



Table 1: Scaling for the synoptic physics described by (4). In the right column the numerical values used. 



2.4 Equations for ^ and r 

Under the boundary conditions (6), the thermal wind relation (^|) takes the simple form = r. 
This identity can be used to eliminate the velocity potential Xi from equations (4.2) and (4.3). This 
results in a closed set of prognostic equations for ^ and r 

d 

—V 2 * = -J(tf, V 2 *) - J(t, V 2 r) - CV 2 (* - r) (7.1) 

d 

— (1 - «V 2 )t = - J(tf, r) + aJ(r, V 2 *) + aJ(#, V 2 t) - h N (r - r*) 

(J L 

- aCV 2 (* - r) + 2aC V 2 r (7.2) 
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where a = a/ fa and r* = O*. In this scaling / 1 « 1.6 • 10 2 is the Rossby number. We will 
study the lowest dimensional nontrivial spectral truncation of these equations. 

2.5 Energy 

In the absence of friction and forcing, the prognostic equations (7) conserve the sum of kinetic and 
available potential energy defined respectively as 

K = a [ f (W • V* + Vr • Vr)dxdy A= f f <d 2 dxdy (8) 
Jo Jo Jo Jo 

in units MD 2 T, 2 . The simplified models will be shown to have a corresponding conserved quantity. 

3 The Galerkin approximation 

In order to approximate equations (7) by a finite number of ODE's we do a Galerkin projection 
onto Fourier modes. On this basis the variables are given by 

^(x, y, t) — ip(m, n, t) exp[i(mfca; + nly)] 

n,va 

t(x, y, t) = ^ 0(m, n, t) exp[i(mfca; + nly)] (9) 

where k = 2ir/s, I — ir. The boundary conditions, and the restriction that 'J and r are real variables 
impose that 

rp(m, —n) = —ip(m, n) 9(m, —n) = —9(m, n) if m =/= 

V>(0,n) = V(0,-n) 0(0, ra) = 6(0, -n) 

ip(m,n) = ip*(—m,—n) 6(m,ri) = 6*(—m,—n) (10) 

This Fourier decomposition is equivalent to the introduction of a basis of eigenfunctions of the 
Laplacian operator on our domain, with the specified boundary conditions. The eigenfunctions are 
4>0n = cosnly for zonal wavenumber zero, and (j> m n — e lmkx smnly otherwise. 

If we apply the zonally symmetric forcing 0* = |AT<^oi, with temperature contrast AT between 
the boundaries, there is an exact solution to equations (7). It is given by 

This solution is called the Hadley state and describes a strong jet in the upper layer, rising air at 
the south boundary and sinking air at the north boundary. The eigenfunctions with nonzero zonal 
wavenumber describe traveling waves which can grow if the Hadley circulation becomes dynamically 
unstable. 

The projection of equations (7) is given by 

Xcdip(c, d) = -\ cd C(ip(c, d) - 0(c, d)) 

+ klJ2 X rs(ps - qr)S(p + r- c)u{q + s - d){^(r, s)^)(p, q) + 0(r, s)6(p, q)} (12.1) 
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X cd 9(c, d) = -h N (9(c, d) - 9*{c, <Q) - aX cd C{^{c, d) - 6{c, d)) + 2aC'X cd 9(c, d) 
+ klJ2 (P s - <l r ) S (P + r ~ c )m(<Z + a - d){l - a([p 2 - r 2 ]k 2 + [q 2 - s 2 ]l 2 )}^{p, q)6(r, s) (12.2) 

pqrs 

where 5 is the Kronecker delta, 9* is the Fourier transform of r* and \i is defined as 



,x(a) = ( e io, "di/ = 
Jo 



1 if a = 

if a is even (13) 
if a is odd 



The eigenvalues of the operators on the left hand side of equations (7) have been denoted by 
X a b = -(a 2 k 2 + b 2 l 2 ) and X ab = 1 - aX ab . 

4 The six dimensional truncation 

The smallest nontrivial truncation of equations (12) has six degrees of freedom. We set s = 2 and 
define 

2=1=2^(0,1) 2/i = 20(O,l) Ti = 20*(O,l) 

x 2 = 2\/2Re^(l, 1) 2/2 = 2\/2Re0(l, 1) T 2 = 2V2Rc(9*(l, 1) 

z 3 = 2\/2Im^(l, 1) 2/3 = 2\/2Im0(l, 1) T 3 = 2\/2Im0*(l, 1) (14) 

These variables satisfy the following equations 



Xi = 


-C(xi - 2/1) 








(15.1) 


X\l±2 = 


-AnC(x 2 - 


2/2) + Ai (5(xix 3 


+ 2/12/3) 




(15.2) 


A11X3 = 


-AhC^s - 


2/3) - Ai <5(xix 2 


+ 2/12/2) 




(15.3) 


A012/1 = 


— aAoi(Ca;i - 


-[C + 2C']»i)- 


h N (yi - 


T\) + S(x 3 y 2 - x 2 y 3 ) 


(15.4) 


A112/2 = 


-aAn(Ca;2 - 


-[C + 2C"]2/ 2 )- 


ftiv(j/2 - 


T 2 ) + <5(Ai xij/3 - v w x 3 yi) 


(15.5) 


AllJ/3 = 


-aAn(Ca;3 - 


-[C + 2C']y 3 )- 


/ijv(2/3 - 


T3) - 5(Ai a;i2/2 - v w x 2 yi) 


(15.6) 



where 8 = 8kl/(3w) and v a \, = 1 + aA a fc. When we put dissipation and forcing to zero, the ODE 
system (15) has a conserved quantity L, defined by 

L = -aXoixj - aAn(a;2 + x\) + A012/1 + An(y| + 2/3) (16) 

which corresponds to the projection of the sum of kinetic and available potential energy defined in 
equation (8). Using L as a Lyapunov function we can show, that a trapping region for equations 
(15) is defined by 

L < 3 " T " 2 (17) 
" {h N - 2aAnC") 2 [ ' 

where ||.|| denotes the L 2 norm, and we have assumed that < Iin < C, < 2C" < C and a < 1, a 
realistic range for the parameters. Note, that the divergence of the vector field, defined by equation 
(15), is constant and negative. Therefore, volume elements always shrink. 
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4.1 Bifurcation diagram 

In figure (||) the partial bifurcation diagram of system (15) is shown. Th is diagram, and all other 
diagrams below, have been calculated using the software package AUTO ( Doedel et al. [1986 1). To 
obtain this picture, we varied T\ and T2, setting T3 = 0. Due to a discrete symmetry in equations 
(15), given by 



x^x' = Rx / 1 \ 

y -> y' = Ry where R= -1 , (18) 

T^T' = RT \ 1 / 

setting T2 = and varying T3 yields the same diagram. This symmetry corresponds to a translation 
x — > x' = x + 1/2 in the physical domain. Furthermore, we have set C = 3.5, = 0.7, C = 0.5 
and a = 0.9. This corresponds to a damping time scale of two days at the earth's surface and 
two weeks at the layer interface. The thermal damping time scale is ten days and the temperature 
difference between the layers is about 34K. 




If we put T2 = and increase T\ from zero, at first a stable equilibrium is the unique limit set 
in the phase space of system (15). This equilibrium corresponds to the Hadley state (|Tl|). At the 
Hopf bifurcation line, marked h, this equilibrium becomes unstable and a stable periodic orbit is 
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created. This is our model's baroclinic instability. The periodic orbit corresponds to a traveling 
baroclinic wave. 

The two line segments fe, joint at cusp point c, denote a fold bifurcation of the equilibrium. 
Within the V-shaped region, bounded by curve fe, there are three equilibria, one of which is stable. 
At the codimension two point, marked fh, the Hopf line and the fold line are tangent. At this 
point an equilibrium exists with one zero e igenvalue and a c omplex pair on the imaginary axis. 



The unfolding of this point can be found in Kuznetsov [1998 1. In the following, we will adopt the 



notation of this book for normal form coeffi cients. An algori thm for computing the normal form 
coefficients of fold- Hop f points is described in Kuznetsov [1999 [ and implemented in the forthcoming 
release of CONTENT ( Kuznetsov and Lcvitin [1997 ]). In this case, we have normal form coefficients 
s = 1 and < 0. 

A torus bifurcation line emerges from point fh, and connects to the flip bifurcation line marked 
fi. Above the torus bifurcation line, the periodic orbit created on Hopf curve h is unstable, below 
the torus bifurcation line it is stable and coexists with a saddle type, two dimensional torus. At 
the point where the torus and the flip bifurcation lines meet, the periodic orbit has two Floquet 
multipliers equal to minus unity. This point is marked ri for 1:2 resonance. At this point, we have 
normal form coefficient s = — 1. Connected to the resonance point ri there is torus bifurcation line 
of the period doubled orbit, which leads to another 1 : 2 resonance point, i'2. In fact, flip bifurcation 
lines fi and f2 are the first two of a period doubling cascade. There seems to be a accumulation 
of 1 : 2 resonance points directly to the right of r 2 . Such an accumulation has been described in 
the context of a biological model in Kuznetsov [1998| , chapter 9.6. Wicczorek et al. [2001] found it 
in a rate equation model for a semiconductor laser. They also provide a partial unfolding of this 
codimension two phenomenon. 

From codimension two point fh two homoclinic bifurcation lines emanate. Both homoclinic 
connections are attached to a saddle focus. Along H 2 the unstable manifold of the saddle focus is 
one dimensional and the saddle value, a, is negative. On this line a stable periodic orbit is created. 
Along Hi, the saddle focus has a complex pair of eigenvalues with positive real part. The saddle 
value is positive along the larger part Hi and another stable periodic orbit is created. However, 
near the leftmost turning point of this curve, the saddle value changes sign. On a small segment it 
is negative, indicating that Shil'nikov type chaos can occur. In the literature, this type of Shil'nikov 
bifurcation, with complex unstable leading eigenvalues, is uncommon. 

The points, where the saddle value is zero, are called neutral saddle focus, or Belyakov, tran- 
sitions. What is known about the unfolding of this transition is summed up in Champneys and 
Kuznetsov [1994 1 . In figure (^) the transition points are shown in more detail, along with a phase 
portrait. An infinite number of cusps of fold lines of periodic orbits are expected to accumulate 
here, corresponding to the creation of successive fol ds of the branch of periodic sol utions which 
becomes homoclinic on Hi. However, as described in Glendinning and Sparrow [1984f[ , the neutral 
saddle focus transition is continuous, and these cusp bifurcations correspond to orbits of very high 
period which are hard to detect numerically. Looking at sections like in figure the case with 
negative saddle values cannot be discerned from the case with positive saddle value. 

In figure (^) a cross section of diagram (0) is shown. We have fixed T\ = and let T 2 vary 
as indicated in diagram (j^) . The top picture shows a number of branches of the period doubling 
cascade, as well as the primary homoclinic branch corresponding to curve Hi. The bottom picture, 
on the same horizontal scale, is a limit point diagram. This picture was obtained by calculating 
the Poincare map on the plane S — {(x, y) 6 K 6 |y 2 = 0}. After a sufficiently long integration, to 
filter out transients, the value of yi was plotted at a number of iterations of the Poincare map. The 
behaviour is chaotic between the accumulation points of the period doubling cascade and changes 
qualitatively near the Shil'nikov bifurcation. 

The lines marked fpj and fp 2 denote fold bifurcations of periodic orbits. Near the cusps of 
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Figure 3: Top: detail of bifurcation diagram (H). The neutral saddle focus transitions have been marked ti t %. To 
the left of these transitions the saddle value is negative and the Shil'nikov condition is satisfied. Bottom: phase 
portrait at (T\,T2) = (0.6,0.3168). These parameter values have been marked with a cross in the top picture. In 
red: the homoclinic orbit Jn blue, dashed: the periodic solution after four period doublings, i.e. the fifth branch in 
the top picture of figure (W). In green: points on the chaotic attractor. 
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Figure 4: Top: continuation of the periodic orbit created on curve h, along the line Tj = TW in fi gure z). 
Solid lines denote stable branches, dotted lines denote unstable branches, Hopf bifurcation point are marked with 
dots. Also shown is the primary homoclinic branch corresponding to curve Hi. Bottom: limit point diagram of the 
Poincare map on S. 
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Figure 5: Cross sections of diagram (Q). Dots denote flip bifurcations, solid lines stable branches and dotted 
lines unstable branches. The labels fi a and fp x 2 refer to diagram (0). Top left: T\ = T< 2 ) = 0.63. The basic 
cycle, born on Hopf line h, and its period double version, are not connected to the homoclinic branch. Top right: 
Ti = T< 3 ) = 0.64. The period doubled cycle now becomes homoclinic on Hi. Bottom left: Ti = T< 4 > = 0.88. The 
branch of the basic cycle folds. Bottom right: Ti = T' 5 ' = 0.98. From one end, Ti increasing from zero, the branch 
of the basic cycle ends in a period halving bifurcation. From the other end, T2 decreasing, the branch of the basic 
cycle becomes homoclinic on Hi. 



these fold lines, the periodic orbit, created at Hopf line h, and its period doubled versions, switch 
branches with periodic orbits that become homoclinic near Hi. This process is illustrated in figure 
(||). In the top left picture, where we have T\ = T^ 2 \ none of the branches of the period doubling 
cascade are connected to the homoclinic branch, see also figure (§)(top). For Ti > however, 
the first period doubled branch becomes homoclinic on curve Hi. The period doubled branch and 
the homoclinic branch collide in a transcritical bifurcation. This process is repeated for the branch 
of the basic cycle, continued from Hopf line h. Therefore, the flip bifurcation lines fi,2,... do not 
simply form a cascade and an inverse cascade. At the fold lines fpj 2 the simple structure of figure 
(|) (top) at Ti = T (1) is rearranged. 

Summarising, the qualitative dynamics of system (15) is as follows. Left of the Hopf line h, 
and within the V-shaped region bounded by fold curve fe, there is a stable equilibrium. To the 
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right of curve h, and below curve fe, the behaviour is periodic before crossing flip bifurcation curve 
fi. When crossing fi near the leftmost fold of the homoclinic curve Hi, where the saddle value is 
positive, a combination of period doubling and Shil'nikov chaos is encountered, as demonstrated 
in figures (||) and ([|). To the right of the neutral saddle focus transition points, the behaviour is 
alternatingly periodic and chaotic. Due to the branch switching, shown in figure (Q), the parameter 
space is divided into small chaotic and periodic windows. 

In figure @ the Kaplan- Yorke dimension of the attractor is shown for Ti = and values of 
Ti at which complex dynamics arise. Most remarkably, the attractor dimension does not exceed 
three. Also, the equilibria and periodic orbits studied in figure ((^) have a feature in common. 
They all have three strongly contracting directions. These observations suggest, that the dynamics 
of system (15) take place on a three dimensional invariant manifold. In the next section we will 
present an approximate, three dimensional invariant manifold of this system, which enables us to 
reduce the model to three degrees of freedom. 



2.4 




2 - 



C 

■2i. 
c 



T3 1.6 



1.4 



1.2 - 



0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 

T 2 



Fi gure 6; Numerical estimate of the Kaplan-Yorke dimension, obtained from an integration during Ai = 1.5 ■ 10^ 
at each parameter value T2, for fixed Ti = T^ 1 ). 



4.2 Reduction to an invariant manifold 



With the parameters set to the values, specified in section (4.1), we can scale the constants in 
equations (15.1)-(15.3) as 



C = eC 



SXiq d 



(19) 
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where e « 1/4 and C w n w 1. System (15) is then written symbolically as 



ex = f(x,y) (20.1) 

y = g(x,y) (20.2) 



where f is defined as 



fx = -C{x x - Vl ) (21.1) 
h = -C{x 2 - 2/2) + k(£i2;3 + yu/3) ( 21 - 2 ) 
h = -C(xa - ya) - k(xxx 2 + 2/13/2) (21-3) 

We assume, that there exists a globally attracting, three dimensional invariant manifold in system 
(20), denoted by W t . This manifolds is represented as the graph of a function, <p e , of the baroclinic 
components: 

W 6 = {(x,y)eIL 6 |x = ^(y)} (22) 

and satisfies 

f(0 £ ,y) = eD0 e - g (0 £ ,y) (23) 



The solution of (|23j) can be approximated asymptotically. Substituting the regular expansion 
4> t — 4>o + etfix + . . . , we find 

1 

Mv)= I A( yi ) B( yi ) |y (24) 



-B{ yi ) A{ yi ) 



where A and B are defined as 



j4(yi) = ^T^f = COS7 and B{vi) = c^y\ =smi - (25) 

This zeroeth order approximation has a clear physical interpretation. It describes a zonally sym- 
metric part of the flow which is equivalent barotropic (i.e. $3 oc ^1) and a phase shift, 7, between 
the traveling waves in the upper and the lower layer. 

The nonlinear terms in equations (15) can be divided into two groups: one represents advection 
of waves with the zonally symmetric flow (the terms proportional to Xx) and the other represents 
energy exchange between waves and the zonally symmetric flow. The nonlinear terms in equation 
(15.4) belong to the second group. If the phase shift is zero, these terms cancel, and the wave 
components are decoupled from the zonally symmetric components. Thus, the function 4>$ describes 
how the energy transfer depends on strength of the zonally symmetric flow. 



The asymptotic expansion of </> e can not accurately describe the solution of equation (23) at 
the realistic parameter value e = 0.25. Therefore, we use a numerical algorithm to continue the 
solution, known analytically in the limit of e 1 0. The algorithm is similar to the graph transform 



described in Brocr ct al. [1997f] . In contrast to the systems considered in their work, however, ours 



is a continuous time system. Therefore, we apply the graph transform to the map, induced by an 
implicit Euler step which approximates the flow of system (20) over a finite time interval. This 



approach was already used by Foias et al. [1988 1 to approximate inertial manifolds by their modified 
Galerkin method. 

Another difference is the choice of coordinates. Here, we do not represent the (approximate) 
invariant manifold as the graph of a section of the normal bundle of some given approximation. 
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Wo 



W:x=<Ky) 



W:x=<Ky) 



Fi gure 7'. Schematic picture of the graph transform. Note, that equilibria of system (20) are intersection points of 
W, W and W e - The dashed lines represent the mesh yijh- 

Instead, we globally represent it as the graph of a function </> £ (y). This can be done provided that 
D0 e has full rank. For small e this condition is satisfied, as we have detD0 £ = 1 + 0(e). When 
increasing e the condition has to be checked numerically. 

Suppose that, for some fixed value of e, we have an approximation, W, of W e . By assumption, 
W e is globally attracting so that the image of W under the flow over a finite time interval of system 
(20) lies closer to W e than W. Thus, by approximating the flow of system (20), we can calculate an 
improved approximation W. For this end we use the map E : K 6 — > M 6 , where (x,y) = £J((x,y)) 
is the solution of 



x = x+— f(x,y) 

e 

y = y + Ag(x, y ) 



(26.1) 
(26.2) 



which defines an implicit Euler time step. The step size, A, is a free parameter. The improved 
approximation is then defined as W = E(W). 

In order to represent W as the graph of a function 0(y), we map a point (<^(y),y) € W onto 
the point (</>(y),y) € W, where <j> is the solution of 



<£ = <£(y- Ag(^,y)) + -f(0,y) 



(27) 



In other words, for a given vector y we look for the point (4>{y'), y') € W which is mapped according 
to (26) onto (<^(y),y) € W. In figure (0) this procedure is sketched. 
Differentiating equation (p7|) with respect to y, we find 



A - - A - 

(I - -D x f + AD0| y ' D K g) T> v 4> = D4>\ y > (I - ADyg) + — D„f 



(28) 
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from which we can calculate T) y 4>, and thus the tangent space y \ W, from T) v 4> and the Jacobian 
of system (20). This enables us to define a local error function. Let F — (f/e, g) denote the vector 
field (20) and V the orthogonal projection onto the tangent space at (</>,y). Then we define the 
error function 

e( y) = ^arccosgglf (29) 
yy) tt |TO,y)|| V ' 

We will represent the approximate invariant manifold on a cubic lattice in IR 3 . The vertices of 
this lattice are located at y^fc = d k) T , where d is the lattice spacing and i, j, k are integers. We 
consider a finite number of lattice points, demanding that all points (4>o{yijk), Yijk) he within the 
trapping region defined in section (Q). In fact, inequality ( |l7| ) is a fairly course estimate. In order 
to reduce computation time and data storage, we use a sharper estimate, obtained by numerical 
computation of the eigenvalues of the linear part of equations (15.1-15.6). To an approximate 
solution <j) we can then assign the error 

£ = max e(yy fc ) (30) 

ijk 

In order to approximate the invariant manifold for finite e, we will proceed in steps. First of 
all, we fix an initial value for e and calculate <pi and D0i on each lattice point. Next, we solve 
equation (p7|), by Newton iteration, and subsequently (p8|), to find the next approximation (j> and 
its derivatives. A suitable initial guess for the Newton iteration is obtained by linearisation in 
£ = y' — y. This yields 

from which we can estimate £ and, subsequently, <j). The graph transform is iterated untill the 
global error is smaller than a fixed threshold £ ma x- If this is achieved, we increase epsilon and 
iterate the process. 

When solving equation (^), evaluation of 4> and D0 inbetween lattice points is necessary. This 
is done by linear interpolation. To this end, each cube in the lattice is divided into six tetraeders 
of equal volume. The associated error is expected to be of order 0(d 2 ). At the edge of the domain 
we consider, it may happen, that evaluation of (j) outside this domain is required. If so, we solve 
equation (|27|) , substituting </>(y') — > 0(y) + D0 • £. The derivatives of 4> are then calculated by finite 
difference. This, however, introduces an error, which disables us to continue the invariant manifold 
up to e = 0.25. 

For example we approximate the manifold at parameter values {T\ 1 T2) — (0.6,0.25). The step 
size is fixed to A ~ 0.0025 and the increment of e is chosen in the range [0.001, 0.005]. The maximal 
error is fixed to £ ma x — 0.05. At each value of e about six iterations of the graph transform are 
needed. The computations were done with a lattice spacing of d = 0.01, in the trapping region 
L < 0.52. Thus, about 5 • 10 5 points on the manifold are calculated. 

In figure (^) the approximate invariant manifold is illustrated. Shown are the stable periodic 
orbit of system (20) at e = 0.1 and (Ti,T2) = (0.6,0.25) and a forward integration of the system 
y = si^eiy), y)- To find the latter integral curve we approximated every point of </> e (y) needed by 
linear interpolation of the results of the calculation presented above. 



5 The reduced system 

In the following, we will describe system (20), reduced to the approximate invariant manifold Wq. 
This reduction can be done analytically, so that we can compare bifurcation diagrams. Physically, 



16 




Figure 8: Time series of the zonally symmetric component, yi, and a wave component, 1/3. Solid: solution of system 
(20) at (T1J2) = (0.6, 0.25) and e = 0.1. Dashed: solution of y = g(</> e (y), y), using the approximate solution <f)c of 
equation (E3) obtained by the graph transform method. 



an argument to study the system reduced to Wq, rather than to the numerically approximated W e , 
is that the quantitative error, introduced by setting e = 0, is smaller than the error introduced 
by the Galerkin approximation. Equations (15) form a qualitative model of one aspect of the 
atmospheric circulation, namely the interaction between the jet stream and the baroclinic waves. 
A further simplification is justified if it keeps the qualitative behaviour in tact. 

Substituting expression (M), we find that the reduced system, y = g(<fio(y), y), is given by 



yi = -c 1 y 1 -d 1 {yl+yl)+T 1 (32.1) 
2/2 = -C22/2 + c 3 y 3 + d 3 yxy 3 + d 2 yiy 2 + T 2 (32.2) 
2)3 = -C22/3 - c 3 y2 - d 3 yxy 2 + d 2 ym + f 3 (32.3) 



where we have introduced 

ci = (~2aX iC + h N )/\ i Ti = h N Ti/X i di = SB/\ 01 

c 2 = {-aXu [C(l -A) + 2C'} + M An f 2 = h N T 2 /\ lx d 2 = 5vB/l n 

c 3 = -aXnCB/Xu f 3 = h N T 3 /Xn d 3 = S[X W - v w A]/Xn (33) 

Again, we have a Lyapunov function 

L = yi + ^(y 2 2 +y 2 3 ) (34) 
^10^01 



and a trapping region defined by 



/ An V IITII 



£ nsscr. |1J - (35) 



In contrast to equations (15), the divergence of the vector field of this reduced system can change 
sign, so that volume elements are not necessarily shrinking. 



As explained in section (4.2), the amplitude of the nonlinear interaction depends on y\ through 
A and B. In figure the dependence of the interaction coefficient in (32.1) is shown. This figure 
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0.2 0.4 0.6 0.8 v 1 1.2 1.4 

y i ^ 

Figure 9: Dashed: the coefficient of interaction, d\, between the jet stream pattern and the wave patterns. Solid: 
the effective dampingcoefficient of 2/2,3 ■ If we fix y\ = y\ = C/k model (32) is equivalent to the Lorenz-84 model 
described in section (0) 



illustrates the life cycle of the baroclinic waves. The linear stability of the Hadley circulation is 
determined by the effective damping coefficient of the wave components 2/2,3, given by C2 — d^yi. 
If this number is positive, the waves are damped. If it is negative, a perturbation of the Hadley 
circulation will lead to growing waves. If y\ is small, there is little energy transfer and the waves 
are damped. As y\ is forced by the meridional temperature gradient, Ti, it grows and the energy 
transfer increases. Beyond y\ ~ 0.276, the effective damping coefficient is negative, and baroclinic 
waves can grow. When y\ = C/k ~ 0.83, the phase shift, 7, is maximal and the waves are optimally 
baroclinic. They extract energy from the jet stream and y\ decreases. Then the waves decay and 
become decreasingly baroclinic in the process. 

In order to see, if system (32) behaves qualitatively the same as system (15), we study the 
bifurcation diagram. It is shown in figure (10). All the bifurcations displayed in diagram (|^) are 
present here, too. Therefore, the qualitative behaviour of the reduced system (32) is the same as 
that of the full system (15). The essence of the extra degrees of freedom in the six dimensional 
model can be captured by the variable coefficients in the three dimensional model. 



6 Reduction to the Lorenz-84 model 

As a final simplification of our model, we fix the coefficients in equations (32). The choice 

A = A(C/k) = B = B(C/k) = 1 

maximizes the efficiency of the energy transfer between the baroclinic waves and the jet stream. 
The phase shift is fixed to 7 = 7r/2. Equivalently, we can reduce model (15) to the tangent space 
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Fi glire 10; Bifurcation diagram of system (32). The enlargement shows the neutral saddle focus transitions ti ; 2. 
Note the similarity to diagram (fcj) 

T( x » jy *)Wo, where (x*,y*) = (C/k, 0, 0, C/k, 0, 0). In other words, we linearise 0o around the 
Hadley state. 

We then scale y and t according to 



"3 



1 r . d 2l c 3 
Vi = -7- C2+C3— \x- — 
d 2 d 3 d 3 



y2 = [c2+c 3 —\—==y y 3 = [c 2 +c 3 — \—^==z 



d 3 ^/did 2 



d 3 ^/d\d 2 



and find 



(36) 



x = —y 2 — z 2 — ax + aF 
y — xy — bxz — y + G 
z = bxy + xz — z 



(37.1) 
(37.2) 
(37.2) 
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Fi gure 1 1 '. Bifurcation diagram of system (37) , with parameters a and b as calculated from the physical parameters 
C, C", /ijv and a. 



which is the model introduced by Lorenz [1984]. For the parameters we find 



a = ci [c 2 + c 3 



d 3 d 3 



G = \J dxd2\C2 + c 3 — ] 2 T 2 

«3 



2 i-i 



(38) 



With the parameters as specified in section (4.1), we thus obtain 



a w 0.35 



1.33 



(39) 



in contrast to the traditional values a — 1/4 and 6 = 4. Diagram ( |ll| ) has been obtained by 
setting the parameters according to (39). There still is a strong similarity to diagrams (Q) and (|l0|). 
However, the neutral saddle focus transitions have disappeared. Along both homoclinic bifurcation 
curves only a stable cycle is created. The accumulation of 1 : 2 resonance points is still there and, 
although homoclinic cur ve H i moves farther away, the branch switching mechanism along fp 2 works 
as described in section (4.1). 
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Figure 12: Bifurcation diagram of system (37), with a = 1/4 and b = 3.24. Torus bifurcation line tr2 and flip 
bifurcation line f2 are tangent at point r. tr2 ends in a 1 : 1 resonance, or Bogdanov-Takens point, marked bt. The 
enlargement shows that extra cusps, C12, have developed on saddle node line fp 2 . Colour coding as in diagram (Till) . 



6.1 Continuation in a and b 

Finally, we continue the bifurcations in diagram ( |ll| ) in parameters a and b in order to establish 
the relation to the diagram at traditional parameter values, presented in [Bhilnikov et al. [1995 ] 



Changing a to its traditional value, a — 1/4, does not change the bifurcation diagram qualitatively. 
When changing b two changes are apparent. Shown in figure dT|) is the dia gram obtained for 
a = 1/4 and b — 3.24. The second torus bifurcation line, tr 2 , connecting two 1:2 resonance points 
in diagram (|Tl|), is now tangent to flip bifurcation line f2 at point r and connects to the fold line 
fp 2 . These two meet in a 1 : 1 resonance points, also called Bogdanov-Takens point. 

Also, a pair of cusp points, marked Ci,2, has developed on fold line fp 2 . These cusps denote cre- 
ation and vanishing of successive wiggles on the branch of periodic solutions connected to homoclinic 
curve Hi. 

Figure (O) shows the bifurcation diagram at parameter values a = 1/4 and 6 = 4. This diagram 



was presented by Bhilnikov et al. [1995]. The torus bifurcation line tr 2 no longer connects to flip 



bifurcation line f 2 and, consequently, has developed an angular degeneracy (Peckham et al. [1995]) 
at point D. Along curve tr 2 the multipliers of the periodic orbit are given by exp(±i(/>), where <p is 
the phase angle. At point D it has a maximum given by <j) 0.87T. 
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Figure 13: Bifurcation diagram of system (37), with a = 1/4 and 6 = 4, see also |3hilnikov et al. 
coding as in diagram (|llh. The angular degeneracy point on tr2 is marked D. 



[1995 . 



Colour 



If the forcing is set to (F,G) — (8,1) the model is known to behave chaotically (e.g. Lorcnz 

[1984| ). An important, unanswered question is how this chaotic behaviour is brought about. With 
the derivation of the model and the bifurcation analysis in mind, we discuss four different possible 
routes to chaos in the following subsections. 

6.1.1 Period doubling cascades 

The flip bifurcation lines / 12 in diagram (|i~3| ) still are the first of an infinite sequence. However, as 
explained in section (4^1), figure (||), a cross section with F fixed does not show a period doubling 



cascade followed by an inverse cascade, due to branch switching along the cusped saddle node 
line fp 2 . The flip bifurcation lines fi,2,..., extending beyond the limits of diagrams (|ll])-(|T^), are 
closed curves. As can be seen in diagram (|ll|), these curves form unnested islands. In accordance 



with the claim to general applicability of the analysis presented in Wieczorek et al. [2001 ], the 
bifurcation diagram of the Lorenz-84 model has much in common with that of their rate equation 
model. Chaotic attractors can be created and destroyed through different routes, including period 
doubling cascades, the Ruclle-Takens scenario and intermittency. 
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6.1.2 Ruelle-Takens scenario 



Another possible route to chaos was proposed by Ruche and Takcns [1971 1. In this scenario, an 
invariant torus is created first. Then a periodic obit appears on the torus, when crossing the 
boundary of an Arnold resonance tongue. If this periodic orbit bifurcates, a chaotic attractor can 
be created. This scenario can be observed in the Lorenz-84 model. We will concentrate on the 
chaotic behaviour at parameter values (F, G) = (8, 1), close to the angular degeneracy point D. 

If we fix F — 8 and increase G from below tr2, we find the attracting period two orbit first, then 
quasiperiodic behaviour in a small interval, then a period 8 orbit on the invariant torus. Figure 
(lif) shows a detail of diagram (|l3|), with the boundaries of the 8 : 1 resonance tongue and two 
subsequent bifurcations, a torus and a flip bifurcation, of the period 8 orbit. 

The phase angle crosses the point <j> — 37r/4 twice, at resonance points Ai,2- One edge of the 
resonance tongue connects these points. It appears that the invariant torus itself goes through two 
subsequent fold bifurcations, so that the edge of the resonance tongue can cross tr2, and in a narrow 
band two stable tori coexist. This is shown in figure (|l5|a). On the inner torus the dynamics is 
quasiperiodic, on the outer torus the period 8 orbits exist. If we further increase G, the period 
8 orbit loses its stability in a torus bifurcation, directly followed by a flip bifurcation, as one of 
the multipliers crosses back into the unit circle through —1. Beyond these bifurcations, a chaotic 
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Figure 15: Poincare sections of the Lorenz-84 model with a = 1/4 and 6 = 4. a: (F, G) = (8,0.9725). Two 
coexisting invariant tori, one with the period 8 orbits, one with quasiperiodic dynamics. The thick dots denote the 
stable orbit, the crosses the saddle type orbit, b: (F,G) = (8,0.99). Beyond the torus and flip bifurcations TR and 
F. Both orbits, marked with crosses and boxes are now of saddle type. 



attractor appears. The corresponding Poincare section is shown in figure (|15|b). The whole picture 
is more involved, as the tori shown in figure (15)a coexist with several tori with phase locked orbits 
of higher period. 

Scenarios for the creation of a chaotic attractor th rough the bifurcation of a periodic orbit on 
on invariant torus are described in Brocr et al. [199c ] . The Ruelle-Ta kens scenario in the vicinity 
of an angular degeneracy was recently found in an electronic model by Algaba et al. [2001 ] . 



6.1.3 Shil'nikov bifurcations 

The neutral saddle focus transitions, found in the six dimensional model, (15), and the approximate 
reduced model, (32), are not present in diagrams (|TT|) - (|T3|) - Therefore, no Shil'nikov type chaos 
occurs in the Lorenz-84 model for these parameter values. We can, ho wever, retrace the tr ansition 
points for different values of b. As the continuation package HomCont ( Doedel et al. [ 1 98^ 1 ) allows 
for three parameter continuation of codimension two points on homoclinic curves, it is possible to 
calculate a curve of neutral saddle focus transitions in the space of parameters 6, F and G. It turns 
out, that two transition points appear on curve Hi if b < 0.419. Just like in diagram (|i~0|), on a 
small segment the saddle value becomes negative. Therefore, for b < 0.419 Shil'nikov type chaos 
can be encountered in the Lorenz-84 model. 



6.1.4 Intermittency 

The last route to chaos described here is through intermittency. If we keep F = 8 fixed and 
increase G, we find a periodic window around G w 1.167. At G = 1.16742, the stable periodic orbit 
u ndergoes a saddle node bifurcat ion. Beyond this point, the behaviour is intermittent. In terms 
of |Pomeau and Manneville 1980 1 , this is type I intermitte ncy. Two time series in the intermittent 



regime are shown in figure (16|). In Wieczorek et al. [2001 1 the same type of intermittency is found 
and pictures of the stable and unstable manifolds of the saddle type orbit are shown. 
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Figure 16: Intermittency in the Lorenz-84 model, a: (F, G) = (8, 1.16743). b: (F, G) = (8, 1.1675). 



7 Conclusion 

The starting point of the analysis in this paper is the truncation to six degrees of freedom of a 
QG two level model of atmospheric flow. Period doubling cascades and ShiPnikov bifurcations 
are identified as routes to chaos in the bifurcation diagram of this model. A measurement of the 
dimension of the chaotic attractor along a section of the bifurcation diagram reveals that it is less 
than three dimensional, hinting at the existence of a three dimensional invariant manifold which 
captures the dynamics. 

In order to approximate this invariant manifold we have introduced a small parameter, e, into the 
equations. In the limit of e J, an exact solution was presented. This solution has a clear physical 
interpretation in terms of energy exchange between the zonally symmetric jet stream mode and 
the traveling waves. Numerical evidence for the existence of this manifold for finite e was found 



implementing a variant of the graph transform algorithm of Brocr ct al. [1997 ]. 

We have shown that we can set e = without destroying the qualitative dynamics. The 
bifurcation diagram of the reduced model in this limit agrees well with the diagram of the six 
dimensional model. The invariant manifold at e = can be linearised around the Hadley state. If 
the six dimensional model is reduced to this linear approximation of the invariant manifold, the 
Lorenz-84 model is found. The parameters of the Lorenz-84 model can then be calculated from 
the physical parameters of the QG model. One of them comes out significantly different. We have 
compared the bifurcation diagram of the Lorenz-84 model at physical parameters to that of the six 
dimensional model and to the diagram at traditional parameter values, first presented in |5hilmkov 
et al. [1995[ | . 

Finally, we have discussed four possible routes to chaos in the Lorenz-84 model. Period dou- 
bling cascades appear both in the six dimensional model and in the Lorenz-84 model. ShiPnikov 
bifurcations are only found in the Lorenz-84 model for parameter values far away from those con- 
sidered here. A route to chaos not considered before in this model is the Ruelle-Takens scenario. 
We have presented evidence that the creation and destruction of invariant tori and the presence 
of resonance tongues lead to chaos in the Lorenz-84 model. Finally, an intermittent transition has 
been presented. The overall picture of chaotic attractors being created and destroyed via different 



routes is reminiscent of the dynamics described in Wicczorck ct al. [2001 1. 

The link between a Galerkin truncation of a QG baroclinic model and the Lorenz-84 model 
justifies the use of the latter in conceptual studies of atmosphere and climate dynamics. It is 
remarkable how much of the bifurcation structure of the six dimensional truncation is preserved, 
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notwithstanding rough approximations (namely setting e = and linearising the three dimensional 
invariant manifold). Probably, this is because the six dimensional model isolates one aspect of 
midlatitude, synoptic flow: the energy exchange between the jet stream and the baroclinic waves. 
When reducing to three degrees of freedom, and subsequently to the Lorenz-84 model, this process 
is modeled qualitatively correctly. 
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